── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.1 ✔ stringr 1.5.2
✔ ggplot2 4.0.0 ✔ tibble 3.3.0
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.1.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)library(plotly)
Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:stats':
filter
The following object is masked from 'package:graphics':
layout
library(knitr)library(zoo)
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Reading and processing the New York Times (NYT) state-level COVID-19 data
1. Read in the data
# load COVID state-level data from NYTcv_states <-read.csv("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv")# load state population datastate_pops <-read.csv("https://raw.githubusercontent.com/COVID19Tracking/associated-data/master/us_census_data/us_census_2018_population_estimates_states.csv")# adjust column namesstate_pops$abb <- state_pops$statestate_pops$state <- state_pops$state_namestate_pops$state_name <-NULL# merge datasets by statecv_states <-merge(cv_states, state_pops, by="state")
2. Look at Data
dim(cv_states)
[1] 58094 9
head(cv_states)
state date fips cases deaths geo_id population pop_density abb
1 Alabama 2023-01-04 1 1587224 21263 1 4887871 96.50939 AL
2 Alabama 2020-04-25 1 6213 213 1 4887871 96.50939 AL
3 Alabama 2023-02-26 1 1638348 21400 1 4887871 96.50939 AL
4 Alabama 2022-12-03 1 1549285 21129 1 4887871 96.50939 AL
5 Alabama 2020-05-06 1 8691 343 1 4887871 96.50939 AL
6 Alabama 2021-04-21 1 524367 10807 1 4887871 96.50939 AL
# format datecv_states$date <-as.Date(cv_states$date, format="%Y-%m-%d")# factor variablesstate_list <-unique(cv_states$state)cv_states$state <-factor(cv_states$state, levels = state_list)abb_list <-unique(cv_states$abb)cv_states$abb <-factor(cv_states$abb, levels = abb_list)# order by state then datecv_states <- cv_states[order(cv_states$state, cv_states$date),]# Confirm the variables are now correctly formattedstr(cv_states)
# Inspect the range values for each variable.head(cv_states)
state date fips cases deaths geo_id population pop_density abb
1029 Alabama 2020-03-13 1 6 0 1 4887871 96.50939 AL
597 Alabama 2020-03-14 1 12 0 1 4887871 96.50939 AL
282 Alabama 2020-03-15 1 23 0 1 4887871 96.50939 AL
12 Alabama 2020-03-16 1 29 0 1 4887871 96.50939 AL
266 Alabama 2020-03-17 1 39 0 1 4887871 96.50939 AL
78 Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
summary(cv_states)
state date fips cases
Washington : 1158 Min. :2020-01-21 Min. : 1.00 Min. : 1
Illinois : 1155 1st Qu.:2020-12-06 1st Qu.:16.00 1st Qu.: 112125
California : 1154 Median :2021-09-11 Median :29.00 Median : 418120
Arizona : 1153 Mean :2021-09-10 Mean :29.78 Mean : 947941
Massachusetts: 1147 3rd Qu.:2022-06-17 3rd Qu.:44.00 3rd Qu.: 1134318
Wisconsin : 1143 Max. :2023-03-23 Max. :72.00 Max. :12169158
(Other) :51184
deaths geo_id population pop_density
Min. : 0 Min. : 1.00 Min. : 577737 Min. : 1.292
1st Qu.: 1598 1st Qu.:16.00 1st Qu.: 1805832 1st Qu.: 43.659
Median : 5901 Median :29.00 Median : 4468402 Median : 107.860
Mean : 12553 Mean :29.78 Mean : 6397965 Mean : 423.031
3rd Qu.: 15952 3rd Qu.:44.00 3rd Qu.: 7535591 3rd Qu.: 229.511
Max. :104277 Max. :72.00 Max. :39557045 Max. :11490.120
NA's :1106
abb
WA : 1158
IL : 1155
CA : 1154
AZ : 1153
MA : 1147
WI : 1143
(Other):51184
min(cv_states$date)
[1] "2020-01-21"
max(cv_states$date)
[1] "2023-03-23"
4. Add new_cases and new_deaths and correct outliers
for (i in1:length(state_list)) { cv_subset <-subset(cv_states, state == state_list[i]) cv_subset <- cv_subset[order(cv_subset$date),] cv_subset$new_cases <- cv_subset$cases[1] cv_subset$new_deaths <- cv_subset$deaths[1]for (j in2:nrow(cv_subset)) { cv_subset$new_cases[j] <- cv_subset$cases[j] - cv_subset$cases[j-1] cv_subset$new_deaths[j] <- cv_subset$deaths[j] - cv_subset$deaths[j-1] } cv_states$new_cases[cv_states$state==state_list[i]] <- cv_subset$new_cases cv_states$new_deaths[cv_states$state==state_list[i]] <- cv_subset$new_deaths}# include in main dataset cv_states$new_cases[cv_states$state==state_list[i]] <- cv_subset$new_cases cv_states$new_deaths[cv_states$state==state_list[i]] <- cv_subset$new_deaths# Focus on recent datescv_states <- cv_states |> dplyr::filter(date >="2021-06-01")# Inspect outliers in new_cases using plotlyp1 <-ggplot(cv_states, aes(x = date, y = new_cases, color = state)) +geom_point(size = .5, alpha =0.5)ggplotly(p1)
p2 <-ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +geom_point(size = .5, alpha =0.5)ggplotly(p2)
# set negative new case or death counts to 0cv_states$new_cases[cv_states$new_cases<0] <-0cv_states$new_deaths[cv_states$new_deaths<0] <-0# Re-calculate `cases` and `deaths` as cumulative sum of updated `new_cases` and `new_deaths`for (i in1:length(state_list)) { cv_subset =subset(cv_states, state == state_list[i])# add starting level for new cases and deaths cv_subset$cases <- cv_subset$cases[1] cv_subset$deaths <- cv_subset$deaths[1]for (j in2:nrow(cv_subset)) { cv_subset$cases[j] <- cv_subset$new_cases[j] + cv_subset$cases[j-1] cv_subset$deaths[j] <- cv_subset$new_deaths[j] + cv_subset$deaths[j-1] } cv_states$cases[cv_states$state==state_list[i]] <- cv_subset$cases cv_states$deaths[cv_states$state==state_list[i]] <- cv_subset$deaths}# Smooth new countscv_states$new_cases <- zoo::rollmean(cv_states$new_cases, k=7, fill=NA, align='right') |>round(digits =0)cv_states$new_deaths <- zoo::rollmean(cv_states$new_deaths, k=7, fill=NA, align='right') |>round(digits =0)# Inspect data again interactivelyp2 <-ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) +geom_line() +geom_point(size = .5, alpha =0.5)ggplotly(p2)
# pop_density vs. casescv_states_today |>plot_ly(x =~pop_density, y =~cases, type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5))
# filter out "District of Columbia"cv_states_today_filter <- cv_states_today |>filter(state!="District of Columbia")# pop_density vs. cases after filteringcv_states_today_filter |>plot_ly(x =~pop_density, y =~cases, type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5))
# pop_density vs. deathsper100kcv_states_today_filter |>plot_ly(x =~pop_density, y =~deathsper100k,type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5))
# Adding hoverinfocv_states_today_filter |>plot_ly(x =~pop_density, y =~deathsper100k,type ='scatter', mode ='markers', color =~state,size =~population, sizes =c(5, 70), marker =list(sizemode='diameter', opacity=0.5),hoverinfo ='text',text =~paste(state, "<br>","Cases per 100k: ", per100k, "<br>","Deaths per 100k: ", deathsper100k)) |>layout(title ="Population-normalized COVID-19 deaths (per 100k) vs. population density for US states",yaxis =list(title ="Deaths per 100k"), xaxis =list(title ="Population Density"),hovermode ="compare")
7. Explore scatterplot trend interactively using ggplotly() and geom_smooth()
Answer: When zooming into September 2021, several Southern states show an increase in naive CFR, reflecting the impact of the Delta wave. The Florida graph shows that the peak in new cases occurred on January 10, 2022, while the peak in new deaths followed approximately 2–3 weeks later. This time delay reflects the progression from infection to severe outcomes and is consistent with expected epidemiological patterns.
# Line chart for naive_CFR for all states over time using `plot_ly()`plot_ly(cv_states, x =~date, y =~naive_CFR, color =~state, type ="scatter", mode ="lines")
# Line chart for Florida showing new_cases and new_deaths together (two lines)cv_states |>filter(state=="Florida") |>plot_ly(x =~date, y =~new_cases, type ="scatter", mode ="lines", name="New Cases") |>add_trace(x =~date, y =~new_deaths, type ="scatter", mode ="lines", name="New Deaths")
9. Heatmaps
Answer: In the heatmap of raw new_cases, large states such as California, Texas, Florida, and New York stand out because their absolute case counts are highest. However, when normalizing by population (newper100k), smaller states with intense outbreaks become more visible. For example, North Dakota, South Dakota, and Louisiana show strong per-capita spikes, while Florida and Mississippi stand out during the Delta surge. This demonstrates how per-capita scaling changes the interpretation of which states were most affected.
# Map state, date, and new_cases to a matrixlibrary(tidyr)cv_states_mat <- cv_states |>select(state, date, new_cases) |>filter(date >as.Date("2021-06-15"))cv_states_mat2 <-as.data.frame(pivot_wider(cv_states_mat, names_from = state, values_from = new_cases))rownames(cv_states_mat2) <- cv_states_mat2$datecv_states_mat2$date <-NULLcv_states_mat2 <-as.matrix(cv_states_mat2)plot_ly(x=colnames(cv_states_mat2), y=rownames(cv_states_mat2),z=~cv_states_mat2, type="heatmap", showscale=TRUE)
Answer: In October 2021, the CFR map shows higher and more uneven fatality rates across states, with certain regions (such as parts of the South (ex. Alaska) and Great Plains (ex. Montana, Wyoming, North Dakota, Idaho)) standing out as darker areas. Within the most recent date, CFR values are lower and more uniform nationwide. This shift reflects improvements in testing, treatment, and vaccination, which reduced both the overall fatality rate and the variability between states.
### For specified datepick.date <-"2021-10-15"# Extract the data for each state by its abbreviationcv_per100 <- cv_states |>filter(date==pick.date) |>select(state, abb, newper100k, cases, deaths) # select datacv_per100$state_name <- cv_per100$statecv_per100$state <- cv_per100$abbcv_per100$abb <-NULL# Create hover textcv_per100$hover <-with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))# Set up mapping detailsset_map_details <-list(scope ='usa',projection =list(type ='albers usa'),showlakes =TRUE,lakecolor =toRGB('white'))# Make sure both maps are on the same color scaleshadeLimit <-125# Create the mapfig <-plot_geo(cv_per100, locationmode ='USA-states') |>add_trace(z =~newper100k, text =~hover, locations =~state,color =~newper100k, colors ='Purples' )fig <- fig |>colorbar(title =paste0("Cases per 100k: ", pick.date), limits =c(0,shadeLimit))fig <- fig |>layout(title =paste('Cases per 100k by State as of ', pick.date, '<br>(Hover for value)'),geo = set_map_details )fig_pick.date <- fig################ Map for today's date# Extract the data for each state by its abbreviationcv_per100 <- cv_states_today |>select(state, abb, newper100k, cases, deaths) # select datacv_per100$state_name <- cv_per100$statecv_per100$state <- cv_per100$abbcv_per100$abb <-NULL# Create hover textcv_per100$hover <-with(cv_per100, paste(state_name, '<br>', "Cases per 100k: ", newper100k, '<br>', "Cases: ", cases, '<br>', "Deaths: ", deaths))# Set up mapping detailsset_map_details <-list(scope ='usa',projection =list(type ='albers usa'),showlakes =TRUE,lakecolor =toRGB('white'))# Create the mapfig <-plot_geo(cv_per100, locationmode ='USA-states') |>add_trace(z =~newper100k, text =~hover, locations =~state,color =~newper100k, colors ='Purples' )fig <- fig |>colorbar(title =paste0("Cases per 100k: ", Sys.Date()), limits =c(0,shadeLimit))fig <- fig |>layout(title =paste('Cases per 100k by State as of', Sys.Date(), '<br>(Hover for value)'),geo = set_map_details )fig_Today <- fig# Plot togethersubplot(fig_pick.date, fig_Today, nrows =2, margin = .05)